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Abstract 

The plasma membrane of human red blood cells consists of a lipid bilayer attached to a regular 
network of underlying cytoskeletal polymers. We model this system at a dynamic coarse-grained 
level, treating the bilayer as an elastic sheet and the cytoskeletal network as a series of phantom 
entropic springs. In contrast to prior simulation efforts, we explicitly account for dynamics of the 
cytoskeletal network, both via motion of the protein anchors that attach the cytoskeleton to the 
bilayer and through breaking and reconnection of individual cytoskeletal filaments. Simulation 
results are explained in the context of a simple mean-field percolation model and comparison is 
made to experimental measurements of red blood cell fluctuation amplitudes. 
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I. INTRODUCTION 

Membranes are essential components of all biological cells In addition to their biolog- 
ical importance, lipid bilayers and biomembranes have also attracted considerable attention 
from physicists due to their fascinating and unusual properties 

as 

. One particularly well 

studied system is the human red blood cell (RBC) membrane. The RBC membrane is a 
composite structure, consisting of a lipid bilayer adhered to an underlying network of fila- 
mentous cytoskeletal proteins (spectrin) via integral membrane protein anchors (see Fig. [1]). 
The spectrin network is observed to be quite regular [4], with an approximate hexagonal 
symmetry extending over the entire cell surface. (It is worth emphasizing that the spectrin 
based cytoskeletal network of RBCs is completely different from the three dimensional actin 
networks common to other types of animal cells 

Given the apparent simplicity of the RBC membrane, it is tempting to attempt modeling 
with elementary elastic models. Indeed, there is a rich history of elastic RBC models to be 
found in the literature 0, E3, 3, ls[ O, [h| . Without attempting a full historical review here, we 
comment that no single elastic model has yet been identified that is capable of reproducing 
the full set of experimental data available for RBC membranes. For example, while the 
work of Lim, Wortis and Mukhopadhyay [8| is capable of capturing the range of observed 
RBC shapes (stomatocyte, discocyte, echinocyte and non-main-sequence shapes as well) seen 
under various chemically induced stresses (e.g. pH, salt, ATP, etc.), this model has not been 
applied to explain the thermal fluctuation amplitudes observed in RBC membranes. And, 
while the models of Gov, Zilman and Safran [7( and Fournier, Lacoste and Raphael |9[ appear 
to do a good job fitting thermal fluctuation data these models do not appear capable 



14j. The failure of 



of explaining mechanical deformation experiments on RBCs 
simple models to consistently explain both thermal fluctuations and mechanical deformation 
experiments has been recognized for years [6|. One recent model does explain both types of 
data within a single elastic model 10|], however this model treats the spectrin network as an 
incompressible and homogeneous viscoelastic plate coupled to the lipid bilayer. It remains 
unclear why such an approximation should suffice for the sparse cytoskeletal network present 
in RBCs. Additionally, this model seems too simplistic to capture the full range of shape 
behaviors explained in reference s|. 

The models mentioned in the preceding paragraph make no mention of the role of energy 
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expenditure in the behavior of RBC membranes. This despite the fact that it is known that 
RBC membranes possess kinase and phosphatase activities capable of altering t he p roper- 
ties of spectrin and other network associated proteins via (de) phosphorylation 15|, [l6j. And, 
certain measurements of RBC fluctuation amplitudes show a correlation between ATP con- 
centration and fluctuation magnitude 0, Q|- One might argue that the difficulty in fitting 
all RBC behavior to a single elastic model stems from the fact that a truly comprehensive 
model must incorporate the effects of energy expenditure by the cell in a realistic fashion. 
Gov and Safran [19( are the first to seriously consider active energy expenditure within the 
RBC from a theoretical standpoint. They have proposed that ATP induced phosphoryla- 
tion and dephosphorylation of the RBC cytoskeletal network leads to a continual dynamic 
evolution of the integrity of the spectrin network. While this picture remains hypothetical, 
without direct proof, it is consistent with the general observations relating ATP concentra- 
tion to membrane fluctuation amplitudes. Gov and Safran (GS) [19j have used this picture 
to motivate a simple picture for RBC fluctuations under the presence of ATP. Local breaking 
and reforming of the spectrin network is captured via non-thermal forces imparted on an 
elastic membrane model. This model has provided the first plausible explanation for the 
viscosity dependence of RBC fluctuations [20]. 

While elastic models with proper accounting for non-thermal energetics may eventually 
prove adequate in describing the long wavelength physics of the RBC, it is clear that wave- 
lengths near or below the spectrin mesh size (~ lOOnm) must be considered within a more 
microscopic picture. A recent model by Dubus and Fournier (DF) 21[ has extended the 
traditional elastic modeling of RBC membranes to explicitly include the cytoskeletal net- 
work at a molecular level of detail. Within this model, the spectrin network is considered 
as a completely regular hexagonal network of phantom entropic springs attached to a fluid 
lipid bilayer. Over wavelengths significantly longer that the spectrin mesh size, the network 
so modeled becomes mathematically equivalent to an imposed surface tension on the fluid 
bilayer. At wavelengths comparable to and shorter than network spacing, the system be- 
haves differently from a simple membrane with applied tension. This model was used to 
compute the spectrum of thermal fluctuation amplitudes for the RBC membrane, but made 
no attempt to account for non-thermal consumption of energy and only computed thermal 
(non-dynamic) observables. 

In this paper, we extend the DF entropic spring model of the cytoskeleton meshwork to 
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include dynamic evolution of the entire system. We allow the anchor points between spectrin 
and membrane to laterally diffuse and we allow for dynamic dissociation and association 
of spectrin links as a molecular level manifestation of the non-thermal energetic picture 
proposed by GS (fig. |2j). 

One important consequence of the GS picture is that sufficiently high ATP concentra- 
tions lead to an appreciable fraction of dissociated spectrin links at the membrane surface. 
Depending upon the timescales for spectrin (re) association, the effective long- wavelength 
elastic properties of the membrane interpolate between two limiting cases. If spectrin 
(re) association kinetics are much slower than all other timescales in the problem, the ef- 
fective tension imposed by the network (as inferred by out-of-plane bilayer undulations) is 
well predicted by a simple percolation-theory argument (see fig. In the opposite limit 
of fast spectrin kinetics, the effective tension is well predicted by assuming each link in the 
network has a reduced spring constant proportional to the probability of the link being in- 
tact at steady state. At intermediate rates, simulations are seen to interpolate between the 
two extreme cases. 

This paper is organized as follows. In Sec. HIl we present our mathematical model for 
the RBC membrane. In Sec. II III details of our simulation methods are discussed. In Sec. IIVI 
results for a fully intact spectrin meshwork are presented, while in Sec. |V]we generalize to the 
more interesting case of a randomly broken network (both static and dynamically broken). 
We discuss our results in relation to experiment in Sec. |VI] and conclude in Sec. IVHI 



II. MODEL 



We treat the RBC membrane as a Helfrich fluid sheet 22j coupled via mobile anchor 
points to a network of springs. Our Hamiltonian is 
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The first term (integral portion) is t 
sheet assuming small deformations 



re standard bilayer bending energy for a Monge gauge 
22| and bending modulus k. A is the projected area 
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of the membrane, r = (x, y) is the position vector in the xy plane and and h(r) is the 
local displacement of the membrane away from the flat reference configuration specified by 
h(r) = (see Fig. [2]). We assume that the lipid bilayer itself has a negligible (bare) surface 
tension. In more general situations, eq. [T] is easily modified to account for non-vanishing 
tension inherent to the lipid bilayer portion of the membrane 0, [3] . The second term (sum 
portion) accounts for the energy of the 2D cytoskeletal meshwork, modelled as a network of 
entropic springs (or in other words, ideal chains of polymers) with effective spring constant 
//(// = 3k B T/£ k £ c with £ k and £ c the Kuhn and contour length of the polymer, respectively). 
In a fully intact cytoskeletal meshwork, all spring end points (nodes) are restricted to lie 
on the surface of the bilayer, so their coordinates are specified by ri and /i(i"i), with indices 
i (or j) labelling different nodes. The sum is over all distinct nearest neighbor node pairs 
(equivalently, over all polymer springs), denoted as The factor is included to 

account for the possibly incomplete connectivity of the network. It is equal to 1 when the 
ij link is connected and otherwise. The constant E reflects all other free energy change 
associated with connecting a detached filament end to a node besides the elastic energy 
of the spring (e.g., the binding energy to the node); we assume this energy is negative 
and significantly larger than thermal energy scales to insure stability of the network in the 
absence of non-thermal energy sources. 

Although our starting point is very similar to the DF model, we emphasize a few key 
differences. In DF, the lateral positions of nodes are held fixed in the geometry of minimum 
energy for a flat bilayer surface; i.e. the r» variables are treated as set constants in DF, 
not as variables capable of influencing the energetics and/or dynamics of the system. This 
approximation renders the Hamiltonian analytically tractable, however it is not immediately 
clear that such a choice fully captures all relevant physics in this system. For example, with 
fixed node positions equally spaced on a regular lattice, the "spring" contribution to eq. 
[H amounts to a finite differenced version of the usual surface tension contribution to the 
Helfrich Hamiltonian. At long wavelengths, eq. Q] with fixed nodes is guaranteed to behave 
as a Helfrich sheet under tension. If the nodes are mobile, as physically expected, the spring 
network represents a simple manifestation a tethered membrane. Such membranes are known 
to exhibit more complicated fluctuations than expected for Helfrich fluid bilayers 23j. We 
also emphasize that much of the following work is concerned with membrane dynamics, 
which was not considered in DF. One of the most interesting aspects of our model is the 
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dynamic breaking and reformation of cytoskeletal filaments, which could not be studied with 
the equilibrium approach adopted by DF. 

Dynamics in our system are overdamped, owing to the low Reynolds number environment 



present at cellular length scales 



241 ] . For bilayer height fluctuations, we have the following 



Langevin type equation of motion, which accounts for hydrodynamic flow in the surrounding 
cytoplasm 25, Q, 3] 

dh(r, t) 



dt 



/oo 
rfr'A(r-r')[F(r',t) + C(r',t)] 
-oo 



Here, A(r) is the diagonal part of the Oseen tensor 28j, given by 

1 



A(r) 



(2) 



(3) 



87r?7|r 

where rj is the viscosity of the surrounding fluid. The above integral is taken over the entire 
x,y plane; it is assumed that the area of interest, A, is embedded within a periodically 
repeating environment of identical subsystems. F(r, t) is the force per unit area on the 
membrane, 

5H 



- J>e y -(^ 2 (r - Ti)[h(T,t) - h( rj ,t)], 

<3> 



(4) 



where 5(r) is the Dirac delta function, and the sum is over all nearest neighbors of node i, 
denoted as (j). ((r,t) is a spatially local Gaussian white noise satisfying the fluctuation- 
dissipation relation 

(CM)> = o, (5) 

(C(r, t)C(r', 0) = 2k B TA-\r - r')6(t - t'), (6) 

where kg is Boltzmann's constant and T temperature of the system. 

We have another set of Langevin equations describing lateral diffusion of the nodes within 
the bilayer 



dt 



D 



(7) 



with 
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^2f4a(t) \(Ti - r,-) + [h(Ti) - h( rj )] 
(j) 1 



dh{ri 
dTi 



(8) 



and 



(CM) = o, ^ (9) 
(Ci(*)-Ci(0> = ^j^SiAt-t'), (io) 

where D is the lateral diffusion constant of the node across the membrane surface. Eqs. [2] and 
[7] completely specify the dynamics of the lipid bilayer and the attached 2D meshwork. Notice 
that the two sets of Langevin equations are coupled (and must be solved simultaneously) 
via the shape of the membrane surface. We note that eq. [8] neglects the purely geometric 



effect of non-flat membrane geometry on the (x, y) motion of node points 29j, |30|]. This 



approximation significantly simplifies our modeling and it has recently been demonstrated 
that such geometric effects are very small for the physical parameters studied herein [scj. 
It is convenient to recast eq. ([2]) in a Fourier basis 2(|. 

^ = A k [F k (t) + &(t)] (11) 

with k = (2TTm/L x ,2irn/Ly) for integer m and n. Here for the latter convenience of the 
simulation, we assume in general a rectangular sample with size L x x L y in real space, and 
therefore the two lattice constants in k space are different. The quantities hk, and 
derive from functions periodic in x and y, due to the assumed periodicity of the system. 
The Fourier transform pair for an arbitrary function, /, with such periodicity is 



= tV 5> eikr ' ( 12 ) 

f k = f drf(r)e-^. (13) 
J A 

The Fourier transformed Oseen interaction, 

f°° 1 

A k = / rfre" ikr A(r) = — , (14) 



in contrast, reflects transformation over the full 2D plane. By construction, the dynamics 
specified by eq. [TT] reflects an infinite network of periodic membrane images interacting via 
the long range 1/r Oseen hydrodynamic kernel. The random forces obey 

(Ck(*)> = 0, (15) 
(Ck(*)Ck'(0> = 2k B TL x L y A k 1 5 K ^6(t~t'). (16) 



For the dynamics of the breaking and reforming of spectrin springs, we consider a simple 
two state kinetic model. We define the rate to reconnect a link as k c , and the rate to 
disconnect a link as kd, irrespective of the instantaneous position of the endpoints of the 
spring. Accordingly, the value of each jumps back and forth between 1 and 0. Defining 
p as the steady-state probability of a link to be connected at any moment, we have 

k c "H kd 

It should be emphasized that our simple model for the breaking and reformation of spec- 
trin filaments does not obey detailed balance since we do not account for the variations 
in energetics caused by positions of the network nodes within the kinetic scheme. This 
approach necessarily corresponds to a non-equilibrium situation, with the dynamics of the 
spectrin network driving membrane fluctuations in a non-thermal manner. Qualitatively, 
this corresponds to the picture proposed by GS, however the detailed kinetics involved in 
the spectrin (re) association process are unknown and likely differ substantially from the 
picture adopted herein. Our simple two-state picture represents one possible manifestation 
of non-equilibrium driving. 

The coupled equations implied by eqs. [11] and [7] are not amenable to analytical solu- 
tion and will be solved via simulation as detailed below. Before proceeding, we note that 
simulations are run on a discrete square lattice. In other words, Eq. f[T5j) is replaced by 

/k = ^a 2 /(r)e- ik ' r , (18) 

r 

where a is the lattice constant and r now only take discrete values on the lattice. Cor- 
respondingly, in k space, the reciprocal lattice (with lattice constant 2n/L x and 2n/L y ) 
contains (L x /a) x (L y /a) points and the summation in eq. [12] is finite. 

Two types of meshworks are considered in our simulations. First, in accordance with 
the true geometry of RBC membranes 3l|, we consider a hexagonal meshwork (Fig. Hk.). 
In such a meshwork, when all links are connected, the average positions of all nodes form 
a hexagonal lattice (this lattice of nodes should not be confused with the lattice used to 
discretize the surface as discussed above). We consider a finite sample and assume periodic 
boundary conditions in a rectangular geometry of size L x x L y approximately commensurate 
with the embedded hexagonal meshwork (see Fig. Hk). The average distance between the 
nearest neighbor nodes, A, is determined from the average surface density of nodes in a 
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RBC membrane. For theoretical interest, we also consider a square meshwork with square 
lattice symmetry for the average positions of nodes (Fig. Hb). Here we simply take a periodic 
square box, i.e., L x = L y = L. 

To connect with experiment and prior theoretical work, our simulations are used primarily 
to calculate the mean square height fluctuation of the membrane surface in k space, (|/ik| 2 ) 
(angular brackets represent non-equilibrium averages as well as equilibrium averages in this 
work). At long wavelengths (A ^> A) we observe that the relation 

= ksTL x L y 
^ k| 1 K eff k* + a eff k2 V ' 

holds fairly well. This expression corresponds to that expected j^, 0] for a thermal fluid 
bilayer sheet with effective bending rigidity K e ff and effective surface tension u e ff-, but we 
stress that its use in interpreting the simulations is empirical. The composite membrane 
surfaces studied in this work can not be regarded as fluid-like due to the assumed connectivity 
of the cytoskeleton matrix. Much of our analysis is presented in terms of & e ff, which is 
obtained by fitting the simulated data to eq. [19] (while assuming n e ff = k unless noted 
otherwise). For future reference, we define the "free fluctuation spectrum" of the sheet, 
(|^k| 2 ) /, as the result anticipated from eq. [T]in the absence of any cytoskeletal contributions 
(i.e. all f y = 0) 

<IM 2 >/ = ^^- (2Q) 

III. SIMULATION METHODS 

Two simulation methods were used to study the composite membrane system: Fourier 



Monte Carlo (FMC) [32l . |33| and Fourier space Brownian dynamics (FSBD) [26l . l34j . 135 ]. 
In the limit of infinitely slow breaking and reformation of spectrin filaments (quenched 
disorder), the ^ variables are static over the course of any finite simulation and the system 
relaxes to a thermal equilibrium dependent upon the connectivity of the network. In this 
limit, both FMC and FSBD simulations may be used to calculate thermal averages and 
both should agree with one another. When spectrin links are breaking and reforming in 
time following the non-equilibrium scheme presented above, we must run dynamic FSBD 
calculations. Our primary interest in this work is the non-equilibrium case, however FMC 
calculations were performed both to verify the accuracy of our FSBD simulations (in the 



9 



static network limit) and to examine the scaling of some our equilibrium results with system 
size (the FMC method is computationally more efficient than FSBD). 



A. Fourier Monte Carlo (FMC) 



The FMC scheme 



32 



is a standard Metropolis algorithm |36j, which uses Fourier 
modes of the bilayer, hk, as the bilayer degrees of freedom. There are two kinds of degrees 
of freedom in our system: membrane shape modes and lateral position of the nodes of the 



spectrin network. For the former, we attempt MC moves on the Fourier modes o: 
to enhance computational efficiency relative to naive real space schemes 
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the system 
. While for 



the latter, we attempt moves that displace the x, y position of the nodes to adjacent sites 
of the real space lattice inherent to the simulations. Note that the shape of the membrane 
surface is continuously variable in our description, while the lateral position of nodes is 
discrete. The primary advantage to evolving the shape of the bilayer surface in Fourier space 
is that we may tune the maximal size of attempted MC moves for each mode separately. In 
the case of a simple fluid bilayer (without attached cytoskeleton) under finite surface tension 
it is clear that a good choice for maximal move sizes is (see Eq. (1) in Ref. 371 ]) 

(nk 4 + ak 2 )5l 



2L x L y 



const., 



(21) 



where 5^ is the maximum attempted jump size of mode k. In our case, a similar choice 
works well only for wavelengths sufficiently long that eq. [19] is obeyed, however it is a simple 
matter to tune each <5k individually to optimize simulation efficiency. In practice, maximal 
jump sizes were tuned to insure that the acceptance ratio of all trials was approximately 
1/2. 



B. Fourier space Brownian dynamics (FSBD) 



26 



351 ] . Here we 



The FSBD method has been fully detailed in our previous work 
present only a minimal discussion to introduce the method. In this section we consider the 
simple case when in Eq. is independent of time and postpone the discussion of Pij(t) 
for the next subsection. 
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Integrating Eqs. (fTTT) and (J7|) from t to t + At for small At, we have 



h k (t + At) = h h (t) + A k [F k (t)At + R^(At)}, 

n(t + At) = n(t) + —\Fi(t)At + Rj(At)]. (22) 
k B T 

where 

/t+Ai pt+At 
dt'C^t'), Ri(At) = y dt'CiC*')- (23) 

The statistical properties of i? k (At) and Rj(At) follow directly those of Ck(^) an d Ci(X)> 

(i2 k (At)> = 0, (24) 

(R k (At)R k ,(At)) = 2k B TL 2 AtA k 5 k ^, (25) 

(Ri(At)> = 0, (26) 

(Ri(At) • Rj(At)) = A(k B T) 2 At/D5 iyj . (27) 

In the simulation, i? k (At) and Rj(At) are drawn from Gaussian distributions with means 
and variances specified above. Since /i k = h*_ k must be satisfied to insure the height field 
of the membrane is real valued, only about half of the Rk need to be generated in each 
time step; the remaining follow via complex conjugation. The precise formulation of this 
statement is somewhat complicated by the finite number of modes in our discrete Fourier 



transform. Readers are referred to ref. [38| for a detailed discussion of how the full set of 
random forces are to be generated while preserving the real valued nature of h(r). 
At each time step of the FSBD simulation, the following calculations are performed: 

(1) Take h(r,t) and r,(t) from the last time step. 

(2) Evaluate F(r,t) and F;(t) using Eqs. gD and (EJ). 

(3) Compute Fk(t) by Fourier transforming F(r,t). 

(4) Draw _R k (At) and Rj(At) from the Gaussian distributions specified above. 

(5) Compute h k (t + At) and r { (t + At) using Eq. (|22D. 

(6) Get h(r, t + At) through inverse Fourier transformation for use in the next iteration. 
There is one complication in step 2 of the above procedure. While we evolve each r^(t) as 

a continuous variable, the height field of the membrane is only specified at points on the real 
space lattice (more precisely, it is only readily obtainable via Fast Fourier Transformation 
at these points). To compute the forces for use in eqs. H] and El both r« and h(ri) are ap- 
proximated by assuming the position of node i directly coincides with the nearest real space 
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lattice site. While this approximation could potentially cause problems due to discontinuous 
jumps in the forces, the surfaces under study are only weakly ruffled. It was verified (in 
the thermal case) that FMC simulations with node positions strictly confined to lattice sites 
and FSBD simulations as outlined here were in good agreement. In practice, we choose At 
small enough that further reduction has no consequence for the reported results, typically 
on the order of 0.01 /is depending upon choice of lattice size a. 



C. Kinetic Monte Carlo (KMC) for spectrin (re)association kinetics 

For the case of a dynamic network, every £ij(t) jumps back and forth between 1 and 0, 
according to the two rates k c and kg (see Sec. [TT]) . We use the stochastic simulation algorithm 
of Gillespie (kinetic Monte Carlo) [39] to pick random times for breaking and reformation 
events to occur. Let us consider the event of reconnecting a link. If the link is disconnected 



at time t = 0, then the waiting time distribution for link reconnection is W(t) = k c e~ kct [40]. 
A random number consistent with this distribution is obtained by applying the following 
transformation to a uniform random deviate x 41| 



t=-—lnx. (28) 

A similar transformation, replacing k c with kd is used to determine the time at which an 
intact filament breaks apart. The set of times so obtained provides a complete trajectory 
for the behavior of each ^ for use in eqs. 0] and [HJ The continuously distributed times are 
rounded off to the nearest time point in the discrete FSBD procedure. We note that our 
model assumes each link must always connect the same two nodes (or be broken) - there is 
no provision for a spectrin filament to dissociate from one node and reconnect elsewhere. 



IV. EFFECTIVE SURFACE TENSION IN THE PRESENCE OF AN INTACT CY- 
TOSKELETAL MESHWORK 

In this section we assume a fully connected cytoskeleton meshwork = 1 for all links 
at all times). In this limit, there are no non-thermal effects and either KMC or FSBD 
simulations may be performed to calculate the equilibrium spectrum, (|/ik| 2 )- The qualitative 
features of this fluctuation spectrum have been predicted by Fournier et. al. [9]. They 
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argued that the behavior of the composite membrane surface should behave simply in two 
limits. In the short wavelength limit, the effect of the cytoskeleton might be expected 
to play a minor role; neglecting the cytoskeletal terms in eq. [T] leads to the prediction 
(|^k| 2 ) = ksTL x Ly/ 'nk 4 . In the long wavelength limit, the cytoskeleton should play an 
important role, but may be regarded as a continuous medium imparting an effective surface 
tension to the bilayer (and possibly modifying the bare bilayer bending rigidity) as in eq. [191 
"Short" and "long" wavelengths referenced above are understood to be interpreted relative 
to the size of the individual spectrin links (A, see figH]). The original work of Fournier et. 
al. assumed a sharp transition between these two regimes and fit experimental data with a 
transition at a wavelength of approximately 2ttA. The simulations below are in qualitative 
agreement with this picture, but place the transition wavelength at A and predict a finite 
width to the crossover region. 

Both FMC and FSBD simulations were performed with identical results (as expected). 
Simulations were seeded from an initially flat membrane with the cytoskeletal anchor points 
arranged in a perfect lattice (as indicated in fig. HJ). The initial configuration was allowed 
to fully equilibrate before collecting any data for analysis. The real-space lattice constant 
used in the simulations was a = A/4 with box dimensions of L x = 32a and L y = 42a for the 
hexagonal network and L x = L y = L = 32a for the square network (except where indicated 
otherwise, these values of a and L x ,L y and L were used in all reported simulations). This 
corresponds to 96 independent nodes (192 triangular corrals within the periodic box) in the 
case of six-fold connected anchors and 64 independent nodes (64 square corrals within the 
periodic box) in the case of the four-fold connected anchors. In preliminary runs, it was 
verified that neither increasing the sample size (L x ,L y ) nor decreasing a by a factors of 2 
significantly altered results; the values outlined above were thus chosen to insure converged 
results with minimal computational expense. For convenience, all physical parameters used 
in the simulations are listed in Table [B 

Our simulation results are shown in Fig. [5j In the long wavelength limit, the results are 
well fit by eq. [[9j assuming K e ff = n and using 

a s = /i, a h = V3fi. (29) 

for the effective tension in the square and hexagonal symmetry simulations, respectively. 
These tension values are not fit constants, but rather may be inferred from the cytoskeletal 
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contribution to eq. [TJ Surface tension is an energy per unit area, so we may calculate its 
value as the ratio of entropic spring (cytoskeleton) energy per corral to the area per corral. 
In the square geometry, there are effectively two springs per corral (each spring is shared by 
two adjacent corrals) and using an idealized zero temperature geometry, a single corral has 
area A 2 and total spring energy of 2 x ^A 2 . The reported value for a s follows immediately. 
A similar calculation leads to the somewhat larger value of cr^ in the hexagonal geometry, 
reflecting the higher density of springs in this case. The numerical value of so calculated 
is 1.3 x 10 3 keT /iin" 2 = 5.3 x 10~ 6 J m" 2 , which is close to a theoretical fit 9j of the 
experimental result [u| (see Fig. [5]). 

The fluctuation spectra in fig. indicate that free membrane predictions hold reasonably 
well out to wavelengths of approximately A (A; -1 = A/(2n) ~ 0.02/xm), with good con- 
vergence to long-wavelength behavior established by 5A (&T 1 ~ O.lyum). The intermediate 
regime between wavelengths of A to 5A encompasses the crossover between the two limiting 
cases. While this fact is unfortunate in light of the experimental data for the RBC 111 ] . 
which displays a crossover at longer wavelengths, the simulation predictions are unambigu- 
ous. The experimental data remains a mystery, but we do note that it may be accounted 
for by adoption of an ad hoc harmonic confining potential Q]. 

One of the motivations for this work was to test the performance of the analytical theory 
developed in DF. In particular, we anticipated that the mobility of cytoskeletal anchors 
would lead to some quantitative deviations from the DF theory at short wavelengths. In 
fact, the mobility of anchor points leads to insignificant changes in the fluctuation spectrum 
for physical parameters relevant to the RBC (see fig. [6]). We also note that the detailed 
analysis of DF does slightly better in reproducing the simulated spectrum than does the 
adoption of surface tensions implied by eq. [29] (see fig. [7j). At intermediate wavelengths, 
the small deviations of K e ff away from k predicted by DF do improve the fit as compared 
to our naive arguments, however the effect is very slight in comparison to the leading order 
effect of introducing a finite surface tension. As a final point, we note that the anisotropic 
nature of the square network over all wavelengths leads to some variance in fluctuations 
for a given magnitude of k depending on the direction taken. In principle, the hexagonal 
network should not suffer from this effect [3l|, however the underlying square lattice taken 
for our simulations does introduce some anisotropy to the spectra; these effects are most 
severe at short wavelengths. The spread of data points plotted in figs. [5] - [7] reflects this 
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directional dependence in the spectra and should not be taken as evidence of statistical 
noise. As previously mentioned, statistical errors are less than the size of the symbols used 
in plotting. 

V. EFFECTIVE SURFACE TENSION IN THE PRESENCE OF A CYTOSKELE- 
TAL MESHWORK WITH DYNAMICALLY EVOLVING CONNECTIVITY 

We now generalize the results of the previous section to include RBC membranes with 
randomly (and dynamically changing) broken cytoskeletal meshworks. As outlined in section 
Hlj we assume the dynamics of spectrin association and disassociation are governed by simple 
rate processes. The rate constants for connecting a broken cytoskeletal filament, k c , and 
breaking an intact filament, kj, are assumed independent of the distance between filament 
end points on the bilayer surface and independent of all other connections within the network. 
This immediately leads to the conclusion that the probability for a filament to be intact at 
a given time is p = k c / (k c + kd)- Equivalently, p is the average percentage of intact filaments 
in the cytoskeletal network. For the moment, we simply take this picture as a hypothetical 
model for non-equilibrium dynamics in the RBC membrane. A discussion regarding the 
connection between this model and experiment will be provided in sec. I VI I 

Our analysis in this section centers around the calculation of cr(p), the effective tension 
of the composite membrane as inferred from long wavelength fluctuations. As indicated by 
the notation, this tension depends on the degree of connectivity within the network. Not 
apparent in the notation is the fact that this tension also depends upon the magnitude of 
the rates k^ and k c (and not just the ratio of them in p) due to the non-equilibrium nature 
of the dynamics. Extraction of a(p) follows the same general prescription as in the previous 
section; the fluctuation spectrum is collected and compared to the empirical result of eq. 
[T9l At the longest wavelengths modeled in the simulations, it is found that eq. [19] does a 
good job of reproducing the simulation data. Obvious theoretical estimates for a e ff = a(p) 
are only available in the limiting cases of very fast spectrin (re) association kinetics and very 
slow kinetics. In general, the effective tension as a function of p (and kinetic rates) must be 
extracted from simulation. 

In the limit that spectrin (re) association rates are much faster than all other time scales in 
the problem, each cyotoskletal link in the network behaves as an intact link with a diminished 
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spring constant. The numerical value for this weakened spring constant is simply the time 
average of the spring as it flips between the two possible values of \i and zero. This picture 
immediately leads to 

< T s{p)=P^i &h(jp) — VSp/j,. (fast spectrin kinetics) (30) 

Similar arguments have been invoked previously to suggest a possible ATP dependence in 



the shear modulus of the RBC membrane 



191 ] . Such an ATP concentration dependence 



may 



provide an explanation for RBC shape changes as a function of ATP concentration [191 ] . 

In the opposite regime, when network connectivity kinetics are far slower than any other 
timescale in the system, the membrane may be regarded as evolving thermally under the 



influence of quenched network disorder. The behavior o 



the system in the quenched disorder 
limit is analogous to the 2D percolation problem [42J. As such, it is expected that the 
effective tension within the system must vanish at some finite critical p = p c . For values of 
p less than p c , no global connectivity within the network exists and, consequently, there is 
no restoring force possible in response to area dilations of the membrane surface. The 2D 
connectivity percolation limit is known to occur at p c = 2/z, where z is the bond valence 
for each node (i.e. z = 4 with p c — 1/2 for the square network and z = 6 with p c — 1/3 for 
the triangular network). Furthermore, within the approximation of the mean field theory, 
it is expected that the decrease in a as p drops from one to p c is linear. That is 

a s (p) = 2^p - 1/2), a h (p) = ^p//(p - 1/3) 

(slow spectrin kinetics) (31) 

for all p > p c and a = for p < p c . A brief justification for the percolation theory results 
quoted above may be found in the Appendix. 

The discussion of "fast" and "slow" kinetics in the previous paragraph was intentionally 
left ambiguous, without specification of what time scales these quantities were to be com- 
pared with. It seems prudent to avoid a detailed discussion of this issue, due to the many 
different dynamic scales in this particular problem, however a crude discussion is appropri- 
ate. Given our focus on long wavelength elastic properties, the most obvious timescale for 
comparison is the membrane relaxation time for the longest wavelength under observation. 



For equilibrium membranes at tension a and bending rigidity k, it is well known 



25fl that 
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relaxation of is exponential with a characteristic time 

^ - Jhk- (32) 

This result follows immediately from eq. [TTJ Considering only the longest wavelength 
modes in our simulations (k = 2n / ' L max with L max = L for the square and L max = L y for 
the hexagonal simulations) and allowing for tensions between zero and the fully connected 
network values, we find that r m falls in the range of 9 x 10 _4 s - 8 x 10 _3 s. While these values 
are only required to hold for homogeneous equilibrium membranes, they were verified to hold 
reasonably well for the non-equilibrium membranes studied here when a(p) (determined from 
the fluctuation spectrum) is naively used in eq. [32j The relaxation rate associated with the 
two-state spectrin dynamics is (k c + k^), which provides a time scale = l/(k c + k^). The 
fast kinetics limit discussed above would be expected to apply for r m ^> r c d and the slow 
kinetics limit would apply for r m <C Tcd- 

FSBD simulations were run including KMC breaking and reforming events in the spec- 
trin meshwork (as detailed in sec. IHIj) . Although no equilibrium can be reached in these 
simulations due to the intrinsically non-equilibrium nature of the simulation, the systems do 
converge to a steady state regime. (|/ik| 2 ) as a function of k were collected, at steady-state, 
for both the square and hexagonal systems described in the previous section. Three differ- 
ent sets of simulations were run in each geometry corresponding to different spectrin kinetic 
regimes. Specifically, kd was chosen to assume the values 10 4 s _1 , 200 s _1 , and 10 s _1 . A 
variety of different connectivity percentages were simulated, spanning < p < 1. Given 
values for kd and p, k c = k^pj (1 —p) follows immediately allowing for complete specification 
of the model. The three choices for kd roughly correspond to the regimes r c d -C r m , r c d ~ r m , 
and r c d ^> r m respectively. 

Figures M and M display our simulation results, plotted in the form of effective surface 
tensions as a function of network connectivity. The tensions are calculated as in the lower 
panel of Fig. [5] using the largest wavelengths available in the simulations. In addition to 
the FSBD/KMC simulations, FMC results are plotted for the quenched disorder case corre- 
sponding to (k c +kd) = 0. In these simulations, an ensemble of random network connectivities 
consistent with the prescribed p values were run without allowing the network connectivity 
to evolve. The results for (|/ik| 2 ) were averaged over the ensemble to obtain cr(p). FMC was 
used in these cases for numerical efficiency and is valid since the evolution under conditions 
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of static network connectivity is purely thermal. The speed of the FMC algorithm allowed 
a set of simulations to be run for larger membrane sizes in order to clarify the role of finite 
size effects. 

The figures clearly display both expected regimes. Fast network reorganization yields 
results in good agreement with eq. [301 whereas slow reorganization yields results consistent 
with eq. [3TJ Intermediate kinetic regimes interpolate between these two limiting cases. 
Finite size effects appear not to play any role, except perhaps for the case of p values very 
close to p c in the quenched disorder case, which is to be expected due to the percolation phase 
transition at this point. At values of intermediate membrane connectivity, it is possible for 
the systems to display a range of effective surface tension values, depending upon the relative 
timescales for spectrin kinetics as compared to membrane relaxation. The implications of 
this fact will be discussed further in the next section. 



VI. CONNECTION TO EXPERIMENTS 

Evidence that RBC membrane shapes and shape fluctuations depend upon non-thermal 
energy expenditure comes from two different experimental sources. Measurement of cell 



shape [43|] and shape fluctuations [17], [18| under a variety of MgATP concentrations suggest 
this fact directly. Increases in MgATP concentration lead to enhanced membrane surface 
fluctuations. Less directly, it has been shown that RBC membrane fluctuation amplitudes in 



the presence of MgATP depend upon the viscosity of the surrounding solvent [20j. Thermal 
properties of a physical system should not be influenced by transport coefficients, such as 
the viscosity, but should depend only upon system energetics via the partition function. 

Biochemically, it is clear that the presence of MgATP leads to the phosphorylation of 
spectrin [151 ] . which has been implicated in the observed shape changes and fluctuation 



amplitudes in RBC membranes under varying MgATP concentrations [43]. Whether or not 
this phosphorylation event (and/or similar events in other molecular components of the 
network) coincides with breaking of spectrin filaments is less clear, however the hypothesis 
that this is in fact the case [l9|] is compelling and served as one of the major motivations for 
this study. In what follows, we assume that phosphorylation of spectrin induces individual 
filaments to dissociate from the network. We stress that this model is hypothetical and, 
indeed, that some studies refute the correlation between spectrin phosphorylation and RBC 
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shape/elastic properties [44 |45|. 

In a naive model for spectrin phosphorylation, we assume the following kinetic equations 
for spectrin dissociation and association with the network as a whole. 

S a + ATP H S d + ADP (33) 
S d + R 2 h Sa + HPO^". 

In the above, S a and Sd stand for the associated and dissociated forms of spectrin, respec- 
tively, and it is assumed that the reactions proceed via catalysis by kinases and phosphatases. 
We assume ATP concentration is held fixed, either by endogenous biochemistry within the 
RBC in vivo or by experimental control in vitro. Using our previously introduced notation, 
we write [S a ] = Cp and [Sd] = C(l—p), where C is the total concentration of spectrin on the 
cell surface. The steady state condition is formulated as dp/dt = 0, which, in conjunction 
with kinetic eqs. [331 leads to 

p = (34) 
F n c + 2n y ' 

Here, both kinetic constants have been wrapped up into the single constant n c . n is the 
concentration of ATP in solution. n c has been defined such that the condition n = n c 
implies p = p c = 1/3 (in this section we assume the hexagonal network due to its direct 
connection to RBC systems) , i.e. n = n c specifies the concentration of ATP necessary 
to reach the percolation phase transition in systems with quenched disorder. While this 
model neglects the observed dependence of spectrin dephosphorylation on n 15| . it has 
the advantage of simplicity in the absence of quantitative kinetic data and follows the line 
of reasoning previously introduced Ji|, allowing for comparison to earlier work. Eq. [3H 
predicts the connectivity of the network as a function of ATP concentration and allows for 
comparison between the modeling of the previous section and experimental data. 

Experimental data for RBC fluctuations as a function of ATP concentration is available 



only in terms of real space height amplitudes, h(r) [l8( (see fig. [10]). Following Eqs. (|T2|) 
and ( TT91) (and for simplicity assuming L x = L y = L) , 

where the interval between adjacent k x or k y is 27r/L. In the limit of large L, the sum 
can be approximated by the integral L 2 /(2ir) 2 J dk, and we have 

(Mr) 2 ) = ^ln(l + ^). (36) 
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The limiting expressions and simulations of the preceding section provide a route toward 
estimation of a(p) under various rates of spectrin association/dissociation. Eq. 1341 enables 
us to translate these results into tensions as a function of ATP concentration. The resulting 
tensions may be used in eq. [33 to calculate the real space membrane fluctuation amplitudes 
as a function of ATP concentration. RBC's have a diameter of 7 microns and we use 



7/iin 



46| in our comparison to experiment. It is also important to note that experimental 



observations of RBC fluctuations were carried out on cells that were allowed to "firmly adhere 
to [a] glass substratum" 18j under the influence of their natural (presumably of electrostatic 



origin) affinity for such surfaces. The adhesion of the bilayer to an underlying supporting 



matrix wi 
RBC (47, 



1 be assumed to contribute a bare surface tension, Ubare to the membrane of the 



481 ] . That is, we assume o = cr bare + cr(p). 



We take n c and <7f, ar . e as two fitting parameters to reproduce the experimental data. The 
resulting fits, assuming the two extreme cases of fast spectrin kinetics (eq. [30] assumed) 
and slow spectrin kinetics (eq. [31] assumed) are displayed in fig. [TO] The values adopted 
by the fitting parameters in these two cases are Obare = 0.2/z = 150kBT//zm 2 (both cases), 
and n c = 0.22 mM for the fast kinetics case and n c = 0.5 mM for the slow kinetics case. 
It turns out that the limiting case of slow spectrin kinetics provides the most satisfactory 
fit of the data, not only in comparison to the fast kinetics case, but also with respect to 
intermediate kinetic regimes; the sharp onset of tensionless fluctuations at n c naturally leads 
to the observed plateau in the experimental data. 

The slowest possible time-scale associated with RBC membrane fluctuations in our model 
(r m ) is on the order of a couple of seconds. This assumes a tension-free membrane and 
L = 7/im. Our data analysis suggests that spectrin breaking/reformation dynamics are 
substantially slower than this, since the quenched-disorder (percolation theory) results pro- 
vide the best fit to the experimental data. This finding would appear to be consistent with 
the experimental observation [ijj] that establishment of steady-state levels of spectrin phos- 
phorylation occurs on the time scale of minutes to tens of minutes with changes in MgATP 
concentration. Previous theoretical work exploring the consequences of spectrin dissociation 



on membrane fluctuations 



191 ] has assumed the time-scale for spectrin kinetics to be on the 



order of milliseconds. Within the context of the present model, rates this fast seem unlikely. 

Experimentally, it is observed that RBC fluctuation amplitudes depend not only upon 
MgATP concentration, but also upon the viscosity of the solvent surrounding the mem- 
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brane [20(| ; increasing viscosity above that of aqueous buffer solution leads to decreases in 
observed fluctuations. Membrane fluctuations at thermal equilibrium should yield identical 
statistics regardless of solvent viscosity (although the dynamics will certainly differ), so these 
experiments provide additional proof of the non-thermal character of RBC fluctuations. 

Within the present model, Eq. (I32p clearly displays the linear relationship between inverse 
viscosity and the relaxation timescales for individual membrane modes. In the preceding 
section, it was argued that the ratio between membrane relaxation rates and spectrin kinetic 
rates determines the magnitude of effective tension that must be used if one desires to 
fit membrane fluctuations to the functional form of eq. [191 At a given connectivity, p, 
fluctuation amplitudes decrease {cr e ff increases) as the spectrin kinetic cycle increases in 
speed. Alternatively, since it is the ratio of kinetics to membrane relaxation that is relevant, 
we expect fluctuation amplitudes to decrease as viscosity is raised assuming all other system 
parameters are held fixed. Qualitatively, this trend is in agreement with the experimental 
results. 

To obtain quantitative comparison with the variable viscosity data it is necessary to 
evaluate (h(r) 2 ) explicitly from direct simulation, without assuming a single effective tension 
as in eq. [35] Within our model, fluctuation amplitudes are affected by viscosity changes only 
by shifting the timescale for membrane dynamics relative to spectrin kinetics. A relative 
shift in timescale (at constant p) corresponds to a vertical motion between the limiting 
regimes plotted in fig. [9j In order for viscosity changes to yield any effect, the timescales for 
spectrin kinetics and membrane relaxation must be comparable, (e.g. if spectrin dissociation 
takes an hour, increasing viscosity by a factor of 10 will not remove you from the limiting 
"slow kinetics regime" . Similarly, if dissociation took a nanosecond, no reasonable viscosity 
change could promote the system out of the "fast kinetics regime".) And, if these timescales 
are comparable, differences in the individual relaxation rates at each wavelength will lead to 
different effective tensions as a function of wavelength. In practice, it is not feasible to run 
simulations for a full 7\im x 7\im patch with sufficient sampling to obtain reliable results. 
It is possible to run a 1/im x 7p,m patch. The effective tensions as a function of wavevector 
were extracted from this rectangular geometry and were used in a generalized version of eq. 
[35] with a — > &(k) to obtain (h(r) 2 ) for the full 7/xm x 7\im system. 

To simultaneously find agreement with the ATP concentration data (fig. [TO]) and the 
variable viscosity data it is necessary to assume spectrin rates that are slow relative to 
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membrane relaxation when the solvent is pure buffer solution (as in the conditions of ref. 
[3] and fig. [TO]) , but are poised to move out of this limiting regime with small viscosity 



increases. By trial and error, it was found that t c ^ = 5 x 10 5 //s for n = 1.2 mM (the ATP 



concentration conditions of ref. 2(|) satisfies this condition and yields the best fit to both 
data sets. Other model parameters were taken identical to the "slow spectrin kinetics" fit 
to the data of fig. [101 In particular, L = 7/im, n c = 0.5 mM and a^are = 150kBT//iin 2 . 
In fig. [TJJwe plot real space membrane fluctuations as a function of rj~ l for viscosity values 
spanning an order of magnitude. We observe a reduction in fluctuation amplitudes of about 
over the entire range of viscosities studied, which falls just outside the errors of the 



experimental data [20|| . It is important to stress that figs. [TD] and [TT] simultaneously fit two 
very different experiments to a single set of physical parameters. Our results agree with 
both experiments quite well. 

VII. DISCUSSION 

From a mathematical standpoint, the results of the preceding section could be viewed 
as a success; we fit the available experimental data with our theoretical model. However, 
the physical implications of the derived fit parameters are troubling. Most worrisome is the 
critical concentration of ATP for onset of the percolation limit, n c = 0.5 mM. This implies 
that at physiological conditions (n ~ 2 mM), the spectrin network is 89% dissociated. A 
similar degree of dissociation (83 %) is predicted for the n = 1.2 mM conditions of the 
variable viscosity experiments described above. As presented, it is clear that our model 
should not be taken seriously in the limit of slow spectrin kinetics and n > n c , since we do 
not allow the connectivity of the network to evolve away from its starting point (i.e. a given 
spectrin tetramer is assumed to always link the same two vertices or be disassociated). In 
a severely compromised network, how could a given filament ever be expected to find its 
assigned association point following a dissociation event? 

As mentioned previously, there does not appear to be definitive experimental data relat- 
ing ATP consumption to _spectrin network dissociation. Rather, this is a plausible hypothesis 



advocated in reference [19|], based upon the limited biochemical data that is available for 
the various constituents of the RBC cytoskeleton. We can not reconcile this picture with 
the simulations carried out in this work, primarily for the reason outlined in the previous 
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paragraph. It is possible that ATP consumption could lead to a significant change in the 
elastic properties of individual spectrin filaments without completely dissociating the fil- 
ament from the network. In such a picture, the two forms of spectrin introduced above, 
Sd and S a , would correspond to filaments with weak and strong effective spring constants, 
respectively (or different natural lengths, etc.). In this context, the analysis we present is 
sensible since network connectivity is always maintained, but the elastic properties of indi- 
vidual links are changing. The percolation limit in this picture corresponds to the point at 
which it is impossible to follow a path across the spectrin network without stepping on at 
least one weak Sd link. Although such a picture is completely hypothetical, we note that 



it is theoretically possible to alter the effective e 



maintaining lateral integrity of the network 
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astic properties of spectrin filaments while 



50J. 



The underlying physical picture pursued in this work was motivated by and is similar in 
spirit to the work of Safran and Gov [3] (i.e. a dynamically breaking and reforming spectrin 
network). However, the simulation model implemented to study this system is more similar 
to the work of Dubus and Fournier (with the additional possibility of dynamically 
breaking bonds). We have interpreted our simulation results within the context of a simple 
percolation theory. While this picture shares some similarity with the analysis of ref. [3] in 
the context of global effects associated with variable ATP concentration, we find differences 
between the present m ode. and r e, Q in the description of .oca. fluctuate dynanhes. 
For example, we find no evidence in our simulations for localized forces of the sort discussed 
in ref. 19J. Instead, the increased amplitudes of fluctuation at high ATP concentration are 
attributable to global properties of the spectrin network within our model. We also find 
that the rates associated with spectrin kinetics must be orders of magnitude slower than 
assumed in ref. [19J in order to reconcile our simulation model with the experiments of refs. 
jisl . Q|. However, it should be pointed out that the particular rates assumed in ref. [3] 
were only rough estimates and are not essential to the qualitative predictions of that work 

H. 

The starting point for all of our simulations, eq. CQ assumes a central-force network for the 
cytoskeletal matrix and a phantom entropic spring description for the spectrin filaments. In 



addition, we do not allow for the possible formation of 5-fold or 7-fold defects [19] , |52| in our 
network. Inclusion of direct interactions between spectrin and the bilayer surface, generaliz- 
ing the description of the network to deviate from the central-force picture, and/or allowing 
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for defects in network connectivity could all potentially alter the quantitative findings of 
this work. The present simulation model was adopted for ease of numerical implementation 
and to facilitate comparison to the existing work of Dubus and Fournier 21[. This study 
provides a detailed analysis of the behavior of a specific model for a composite membrane 
system (lipid bilayer plus actively labile cytoskeletal network). This model is motivated by 
our (incomplete) picture for the structure and kinetics of the RBC membrane surface and 
may prove useful in the development of more refined models in the future. In particular, it 
should be noted that more refined descriptions of the spectrin network have been developed 
that allow for interactions between the bilayer surface and the filaments and a more complete 



description of polymer elasticity beyond the simple Gaussian-chain model [53j, [5J, [55|. Fu- 
ture modeling, incorporating some of these approaches for the behavior of the cytoskeleton, 
would be highly desirable and is currently under investigation. 

We note in closing that a very recent study by Evans et. al. 56J directly contradicts the 
experimental results in refs. fla. I^. This recent series of experiments finds no correlation 
between ATP concentration and RBC membrane fluctuations. Given the severe disagree- 
ment between experimental results on the RBC system, one has almost complete freedom 
in interpretation of our simulation results. The picture outlined above suggests that we can 
obtain quite good agreement between simulation and experiments that do measure ATP 
dependence of the RBC fluctuations. Within this picture, it seems necessary to assume that 
our "broken" spectrin links actually correspond to intact links with diminished spring con- 
stant and/or a non-zero natural length. Another, equally valid, interpretation of our results 
is that it is impossible to reconcile an actively breaking cytoskeletal meshwork picture with 
available experimental data due to the high degree of network dissociation predicted under 
physiological conditions within such a picture. Given the uncertain experimental landscape, 
it seems impossible to make a definitive statement in favor of either viewpoint. 
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Appendix 



The calculation of an effective tension in the limit of quenched disorder within the spectrin 



network is closely related to the 2D bond percolation problem 



571 ] . A mean field treatment 



for the percolation problem was first developed in the context of electric conductivity 42] 



and was 
finite 
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ater extended to calculate the elastic moduli of networks of Hookean springs with 



60, 



6JJ. The below discussion summarizes the 



591 ] and zero natural lengths 
arguments of Ref. [5^] for the readers' convenience. 

Let us consider a randomly connected meshwork with only a fraction (p) of the links 
intact. Each intact link is a spring with spring constant /i and a natural extension of 
zero. In the spirit of mean field theory, we assume that the global elastic properties of this 
imperfect meshwork are well approximated by a complete meshwork comprised of springs 
with spring constant \i m at every link (see Fig. [T2]) even though individual links will be 
either completely intact or fully broken. Our problem is to determine an expression for fi m . 
For this purpose, let us consider an arbitrary link AB in the meshwork. The real spring 
constant, fiAB, associated with this link is either \x if the link is connected, or if it is not. 
If we apply a force f m across this link (see Fig. [T2"]) . the extension xab is given by 

f„ 



fJ'AB + Heff 



(37) 



where /i e // is an effective spring constant reflecting the contribution of all network compo- 
nents except the link AB. It can be shown that 



58] 



fji e ff = (z/2 - l)/i r 



(38) 



where z is the connectivity of the meshwork, for example, z = 4 for the square meshwork 
and z = 6 for the hexagonal meshwork considered in this work. 

On the other hand, if we assume the mean spring constant \i m for link AB, we predict a 
mean extension, 

x m = • (39) 



To achieve consistency within the mean field approximation, we must have 



Pfm 
H + fleff 



(xab) 

[} -P)fm 

Mf 
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X;- 



(40) 



fn 



+ Heff 



where the second equality is simply a more explicit version of the first. This equation is 
readily solved to yield 

= Z. //, (41) 

1 ~Pc 

where 

Pc = 2/z. (42) 

To calculate cr(p), we simply replace /x in eq. [SUlby fi m from eq. 2U since the latter is now 
the average spring constant of the equivalent complete meshwork. Taking z = 4 and z = 6 
for the square and hexagonal meshwork respectively, we arrive at eq. [311 

One should note that eqs. [31] and 02] are only results of a mean field approximation. In 
the vicinity of p c , it is well known that the correlation effects are strong and that mean field 
theory breaks down. However, the mean field results are adequate to describe the behavior 
of o outside the immediate vicinity of p c as evidenced by figs. M and EJ which is sufficient 
for the purposes of this work. 



[1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of 

the Cell (Galland, New York, 2002). 
[2] S. A. Safran, Statistical Thermodynamics of Surfaces, Interfaces, and Memebranes (Westview 

Press, Boulder, CO, 2003). 
[3] Structure and Dynamics of Membranes, edited by R. Lipowsky and E. Sackmann (Elsevier 

Science, Amsterdam, 1995). 
[4] S. Liu, L. Derick, and J. Palek, J. Cell. Biol. 104, 527 (1987). 
[5] F. Brochard and J. F. Lennon, J. Phys. (Paris) 36, 1035 (1975). 
[6] M. A. Peterson, Phys. Rev. A 45, 4116 (1992). 

[7] N. Gov, A. G. Zilman, and S. Safran, Phys. Rev. Lett. 90, 228101 (2003). 
[8] G. H. W. Lim, M. Wortis, and R. Mukhopadhyay, Proc. Nat. Acad. Sci. USA 99, 16766 
(2002). 

[9] J.-B. Fournier, D. Lacoste, and E. Raphael, Phys. Rev. Lett. 92, 018102 (2004). 
[10] S. B. Rochal and V. L. Lorman, Phys. Rev. Lett. 96, 248102 (2006). 
[11] A. Zilker, H. Engelhardt, and E. Sackmann, J. Phys. (Fr.) 48, 2139 (1987). 
[12] D. Discher, N. Mohandas, and E. A. Evans, Science 266, 1032 (1994). 

26 



[13] V. Heinrich, K. Ritchie, N. Mohandas, and E. Evans, Biophys. J. 81, 1452 (2001). 

[14] H. Engelhardt, H. Gaub, and E. Sackmann, Nature (London) 307, 378 (1984). 

[15] W. Birchmeier and S. J. Singer, J. Cell Biol. 73, 647 (1977). 

[16] V. Bennett, Biochim. Biophys. Acta 988, 107 (1989). 

[17] S. Levin and R. Korenstein, Biophys. J. 60, 733 (1991). 

[18] S. Tuvia, S. Levin, A. Bitler, and R. Korenstein, J. Cell Biol. 141, 1551 (1998). 

[19] N. Gov and S. A. Safran, Biophys. J. 88, 1859 (2005). 

[20] S. Tuvia, A. Almagor, A. Bitler, S. Levin, R. Korenstein, and S. Yedgar, Proc. Natl. Acad. 

Sci. USA 94, 5045 (1997). 

[21] C. Dubus and J.-B. Fournier, Europhys. Lett. 75, 181 (2006). 

[22] W. Helfrich, Z. Naturforsch. pp. 693-703 (1973). 

[23] D. R. Nelson and L. Peliti, J. Phys. (Paris) 48, 1085 (1987). 

[24] E. M. Purcell, American Journal of Physics 45, 3 (1977). 

[25] R. Granek, J. Phys. II (Paris) 7, 1761 (1997). 

[26] L. C.-L. Lin and F. L. H. Brown, Phys. Rev. E 72, 011910 (2005). 

[27] F. L. H. Brown, Annual Review of Physical Chemistry 59, 685 (2008). 

[28] M. Doi and S. Edwards, The Theory of Polymer Dynamics (Clarendon, Oxford, 1986). 

[29] E. Reister and U. Seifert, Europhys. Lett. 71, 859 (2005). 

[30] A. Naji and F. L. H. Brown, J. Chem. Phys. 126, 235103 (2007). 

[31] D. H. Boal, Mechanics of the cell (Cambridge University Press, Cambridge, UK, 2002). 

[32] N. Gouliaev and J. F. Nagle, Phys. Rev. E. 58, 881 (1998). 

[33] N. Gouliaev and J. F. Nagle, Phys. Rev. Lett 81, 2610 (1998). 

[34] L. C.-L. Lin and F. L. H. Brown, Biophys. J. 86, 764 (2004). 

[35] L. C.-L. Lin and F. L. H. Brown, Phys. Rev. Lett. 93, 256001 (2004). 

[36] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, San Diego, 

2002), 2nd ed., pp. 525-532. 

[37] D. Bouzida, S. Kumar, and R. H. Swendsen, Phys. Rev. A 45, 8894 (1992). 

[38] L. C.-L. Lin and F. L. H. Brown, J. Chem. Theory Comp. 2, 472 (2006). 

[39] D. T. Gillespie, J. Comput. Phys. 22, 403 (1976). 

[40] N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amster- 
dam, 1992), pp. 63,83,220-221. 

27 



[41] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C 

(Cambridge University Press, Cambridge, 1994). 
[42] S. Kirkpatrick, Rev. of Mod. Phys. 45, 574 (1973), and references therein. 
[43] M. P. Sheetz and S. J. Singer, J. Cell Biol. 73, 638 (1977). 
[44] V. P. Patel and G. Fairbanks, J. Cell. Biol. 88, 430 (1981). 

[45] L. Bordin, G. Clari, I. Moro, F. D. Vecchia, and V. Moret, Biochem. Biopbys. Res. Comm. 
213, 249 (1995). 

[46] H. Lodish, D. Baltimore, A. Berk, S. L. Zipursky, P. Matsudaira, and J. Darnell, Molecular 

Cell Biology (Scientific American Books, New York, 1995), 3rd ed. 
[47] F. Brochard, P. G. DeGennes, and P. Pfeuty, J. Phys. (Paris) 37, 1099 (1976). 
[48] L. C.-L. Lin, J. T. Groves, and F. L. H. Brown, Biophys. J. 91, 3600 (2006). 
[49] N. S. Gov, New Journal of Physics 9, 429 (2007). 
[50] Q. Zhu and R. J. Asaro, Biophys. J. 94, 2529 (2008). 
[51] N. S. Gov, personal communication. 

[52] H. S. Seung and D. R. Nelson, Phys. Rev. A 38, 1005 (1988). 
[53] D. H. Boal, Biophys. J. 67, 521 (1994). 

[54] J. Li, G. Lykotrafitis, M. Dao, and S. Suresh, Proc. Nat. Acad. Sci. USA 104, 4937 (2007). 
[55] Q. Zhu, C. Vera, R. J. Asaro, P. Sche, and L. A. Sung, Biophys. J. 93, 386 (2007). 
[56] J. Evans, W. Gratzer, N. Mohandas, K. Parker, and J. Sleep, Biophys. J (2008), in press. 
[57] D. Stauffer and A. Aharony (1994). 

[58] S. Feng, M. F. Thorpe, and E. Garboczi, Phys. Rev. B 31, 276 (1985). 

[59] S. Feng and P. N. Sen, Phys. Rev. Lett. 52, 216 (1984). 

[60] W. Tang and M. F. Thorpe, Phys. Rev. B 36, 3798 (1987). 

[61] W. Tang and M. F. Thorpe, Phys. Rev. B 37, 5539 (1988). 

[62] M. Tomishige, Y. Sako, and A. Kusumi, J. Cell Biol. 142, 989 (1998). 



28 



TABLE I: Physical parameters in simulations 



Parameter 


Description 


Value 


Reference 


K 


bending modulus 


4.8 k B T a 


[5] 


V 


cytoplasm viscosity 


1.5 IcbT s /im" 3 


[5] 


D 


diffusion constant of nodes 


0.53 fim 2 s- 1 


[62] 


A 


the average distance between two nearest nodes (see Fig. [4]) 


0.112 /im 


[31] 




contour length of the spectrin tetramer 


0.2 /im 


[31] 


4 


Kuhn length of the spectrin tetramer 


0.02 /im 


[3JJ 




spring constant of the spectrin tetramer 


750 k B T /im~ 2 


6 



a T= 300K. 

b Calculated using the Gaussian chain model through two parameters above. 
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FIG. 1: A schematic rendering of the RBC membrane, (a): Top View: Looking down on the 
cell with the lipid bilayer invisible for clarity. The nodes of the cytoskeletal network are protein 
complexes embedded in the lipid bilayer. The links are spectrin tetramers, with approximate 
lateral end-to-end distance of lOOnm between the nodes, (b): Side View: The protein anchors are 
confined to the lipid bilayer, but may diffuse laterally. The coils underneath the bilayer are spectrin 
tetramers. The contour length of the tetramers is approximately 200nm, so there is considerable 
extra length coiled up below the membrane surface, (c): One end of a link is disconnected from 
the node, possibly with the help of ATP. 
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FIG. 2: (Color online). A snapshot the RBC membrane as treated in our simulations. The lipid 
bilayer is connected to spectrin via a series of laterally mobile anchoring proteins (white dots). The 
anchors are interconnected via harmonic potentials (not shown for clarity), modeling the entire 
cytoskeletal network as a series of interconnected entropic springs. Our model does not account for 
any steric repulsion between cytoskeleton and bilayer surface, nor does it account for inextensibility 
of spectrin beyond its natural contour length; the individual spectrin links are treated as simple 
Gaussian chains. Individual links in the network are allowed to break and reform following kinetic 
equations that do not obey detailed balance; this provides a model for energy (ATP) consumption 
at the membrane surface. A is the average mesh size (~ lOOnm). 
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(a) (b) 



FIG. 3: A schematic illustration of a randomly broken spectrin meshwork. The solid lines are 
connected links while the dotted lines represent disconnected ones. The percentage of intact links 
(equivalently, the probability for any single link to be intact) will be denoted "p" . All lines together 
(solid and dotted) correspond to a complete hexagonal meshwork. (a): Above the percolation 
threshold, p > p c , the connected links have global connectivity, i.e., they form a cluster extending 
across the entire cell. The corresponding effective surface tension will always be finite, (b): Below 
the percolation threshold, p < p c , there is no global connectivity in the meshwork, i.e. connected 
links form finite islands of connectivity that do not span the cell. The corresponding effective surface 
tension becomes zero if the links are breaking and reforming slowly, (c): Even at instantaneous 
connectivities below the percolation threshold, there can be a finite surface tension if spectrin 
breaking/formation kinetics are sufficiently fast. In this limit, each link can be thought of as a 
spring with a reduced restoring force since out of plane membrane undulations evolve more slowly 
than the lifetime of the link. 




FIG. 4: A schematic illustration of the connectivity of the nodes of the meshwork. (a): the 
hexagonal meshwork. The dotted lines are a guide of eyes for the rectangular sample used, (b): 
the square meshwork. 
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FIG. 5: UPPER: Spectrum of fluctuations, (|/ik| 2 )> (normalized by (|/ik| 2 )/) versus inverse 
wavenumber for the our RBC membrane model assuming a perfectly intact cytoskeletal network. 
The squares and triangles indicate simulation results for square and hexagonal meshworks respec- 
tively. The error bars are approximately the same size as the symbols. The solid lines are theoretical 
results using Eq. (j29|) for a e fj and K e ff = n in Eq. [191 The wavevector k m = 2ir/A is indicated 
on the horizontal axis to show where crossover from free membrane to long- wavelength behavior 



1J. LOWER: 



might be expected to occur. The circles are experimental results taken from Ref. 
The approach of a to its limiting long-wavelength value is displayed. In this case, each point in 
the spectrum on the left was individually fit to eq. [19] by adjusting c e // (assuming K e ff = K )- The 
resulting cr e ff as a function of k~ l are shown. At the longest wavelengths simulated, the inferred 
tensions converge to the predicted values. 
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FIG. 6: Simulation results of (|/ik| 2 ) f° r a hexagonal lattice. The circles are the data for the case 
of immobile nodes placed on a perfect lattice. While the triangles are the data for the case of 
mobile nodes (identical data to fig. [5]). Mobility of the cytoskeleton attachment points does not 
significantly alter the simulation results. 
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FIG. 7: Comparison of different theories for a hexagonal lattice. The circles are the simulation 
results. The dashed line is the prediction introduced in fig. [5j The solid line is the DF model 
prediction. At intermediate wavelengths DF does a slightly superior job in fitting the simula- 
tions. Both approximations fail at small wavelengths and converge to the same behavior at long 
wavelengths. 
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FIG. 8: <j(p)/cr(p = 1) versus p is plotted for a square meshwork. The circles, triangles, and 
diamonds are FSBD/KMC results for kd = 10 4 s~ 1 (r^ <C r m ), kd = 200s" 1 (r c d ~ r m ) and 
kj, = 10s -1 (r c d 3> T m ) respectively. The squares and the stars are FMC results (kd = k c = 0). All 
simulations use box sizes identical to those introduced in the previous section, except for the stars. 
Stars represent FMC simulations on systems that are twice as large in both linear dimensions. For 
clarity, only error bars for the circles and the stars are shown, but all simulations have similar 
errors. The dashed line is the prediction of Eq. ([30]) . The solid line is the prediction of Eq. (f3Tj) . 
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FIG. 9: a(p)/a(p = 1) versus p is plotted for a hexagonal meshwork. All symbols and lines have 
the same meaning as in Fig. only the meshwork connectivity has been changed. (The obvious 
discrepancy between simulation and theory at p = 1 is due to the rectangular geometry of our 
simulation box. It is impossible to perfectly fit a portion of a hexagonal network within a rectangle, 
so the theoretical results that assume perfect hexagonal symmetry are only approximately valid in 
this case.) 




FIG. 10: Real space RBC height fluctuations, \J (h(r) 2 ), as a function of ATP concentration 
(normalized by the zero ATP values, yj (h(r) 2 ) Q ). The circles with error bars are the experimental 
data [18]. The solid line is a theoretical fit assuming ah in Eq. (|3ip . corresponding to the case of 
T m <C T c d (slow spectrin kinetics). The dashed line is a theoretical fit assuming ah in Eq. ([30]) . 
corresponding to the case of t C( i <C r m (fast spectrin kinetics). 
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FIG. 11: Normalized height fluctuations, \J (h(r) 2 ), as a function of inverse viscosity of the solvent. 
The circles with error bars are the simulation results assuming L = 7 fim, cj}j are = 150kBT//im 2 , 
n c = 0.5 mM, n = 1.2 mM and t cc i = 5 x 10 5 /is. The solid line is a fit to our simulation data. 
The dashed line and the two dashdotted lines are the experimental fitting line and the range of 
experimental errors taken from Fig. 2 of Ref. [20I . 




FIG. 12: (a) A randomly broken network of springs is approximated by a complete meshwork of 
springs with diminished spring constant jx m . Here we imagine an attempt to extend the link AB 
using a force of f m . (b) The force f m is resisted both directly by link AB and the rest of the 
meshwork. The latter can be considered as a single spring with effective spring constant /x e //, 



which acts in parallel with AB. (Adapted from ref. 
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